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We show that multiple filamentation patterns in high-power laser beams, can be described by 
means of two statistical physics concepts, namely self-similarity of the patterns over two nested 
scales, and nearest-neighbor interactions of classical rotators. The resulting lattice spin model 
perfectly reproduces the evolution of intense laser pulses as simulated by the Non-Linear Schrodinger 
Equation, shedding a new light on multiple hlamentation. As a side benefit, this approach drastically 
reduces the computing time by two orders of magnitude as compared to the standard simulation 
methods of laser hlamentation. 


The non-linear Schrodinger equation (NLSE), origi¬ 
nally emanating from quantum mechanics, is paradig¬ 
matic of a universal equation which is widely used in 
a variety of fields such as non-linear optics [l|, Bose- 
Einstein condensates Hi, plasma physics 0, or fluid 
mechanics Its analytical properties are quite well- 
known, and exhibit features such as integrability in one 
dimension or finite-time blow-up for higher spatial 
dimensions 0 . 

In the field of non-linear optics, the NLSE describes 
light filaments [1113 forming in the propagation of laser 
pulses which power exceeds a certain critical value. For 
powers much beyond the latter, the beam breaks up 
into many cells, each generating one filament 111- Ijl. 
forming complex multiple filamentation patterns 14| . 
We recently showed that the formation of such patterns 
from an initially smooth laser beam profile defines a 
two-dimensional phase transition governing the geomet¬ 
rical structuring of the beam and the self-organization 
of light filaments 1^. The patterns associated to this 
phase transition are similar to those produced by perco¬ 
lation l3, n I or spin models from the statistical physics 
literature 

The salient features of such systems generally arise 
from the nearest-neighbor interactions between the un¬ 
derlying constituents, mainly quantum or classical spins. 
However, the description of multiple filamentation pat¬ 
terns as the result of basic interacting elements like spins 
was never considered until now. Laser filaments have 
been shown to laterally interact with their neighbors lo¬ 
cated at a distance of several millimeters in the beam 
profile 21 - 23 . This interaction is attractive if the fil¬ 
aments are in phase, and repulsive if they are in an¬ 
tiphase 27, [ill, because it is mediated by interference 
of the photon bath surrounding each filament [29l - l32j . 
However, such interactions have up to now been only 
considered locally. No impact on the global beam profile 
evolution was investigated, or even expected. 

In this Letter, we derive a model for laser multiple 
filamentation, showing that this physical phenomenon 
can be understood as a consequence of self-similarity and 
nearest-neighbor interaction between coarse-grained light 


elements. This results in a description highly reminiscent 
of the Edwards-Anderson spin-glass model , quan¬ 

titatively bridging non-linear optics to out-of-equilibrium 
statistical physics. 

In the following, we will first discuss the self-similarity 
of multiple filamentation patterns. Then, we will show 
that it allows to drastically coarse-grain the dynamics 
with minimal loss of information, provided time is ade¬ 
quately rescaled to account for the change in the speed 
of transverse information flow induced by this procedure. 
The resulting lattice spin model will then be validated 
by a direct confrontation to the results of the standard 
NLSE integration, showing an amazing agreement. 

The starting point of our derivation is the NLSE, 
which, in dimensionless units, reads 


iS^V’ + A'0 + = 0, 


( 1 ) 


where r] is the propagation distance, A = -\- dy the 

two-dimensional transverse Laplacian which accounts for 
geometrical diffraction, and the function / describes the 
nonlinear physical mechanisms at play, including dissi¬ 
pation and saturation. Although the NLSE is ubiquitu- 
ous in physics, in the following we mainly focus on the 
case of multiple filamentation, where pattern formation 
is best characterized both experimentally 37| and the¬ 
oretically [i3|. In filamentation, ip is the electric field 
envelope, and for numerical simulations, it is quite com¬ 
mon to model the non-linearity as 




a2K-2 


( 2 ) 


where the first term accounts for the Kerr self-focusing 
effect, and the two last ones model defocusing by free 
electrons as well as losses due to the K-photon ioniza¬ 
tion releasing these electrons. Without these last two 
terms, some initial conditions of Eq. (HD exhibit finite¬ 
time divergence [s^ - 

Equation HD features a linear instability, called the 
modulational instability, with spectacular experimental 
consequences ranging from the emergence of solitons in 
Bose-Einstein Condensates [s^, 113 to the formation of 
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multiple filamentation patterns ill- 44 1 in large high- 
power laser beams. The growth rate 7 of this insta¬ 
bility can be obtained analytically. For a plane wave 
steady-state writing kj_ the spatial transverse 

wave-vector of the perturbation leads to 221 


7 = fcuY^2V’g/'(V’o) - (3) 

Figure [1] displays the resulting patterns in the case of 
laser propagation in air by solving numerically Eq. O- 
The initial condition is taken as a fourth-order super- 
Gaussian of 5 cm diameter, holding 50 TW at a wave¬ 
length of 800 nm. The relationship between the dimen¬ 
sionless units and the real physical parameters is given in 
the supplementary information [45l |. The modulational 
instability, seeded by the initial beam noise (Fig. [T^), 
triggers the emergence of mesoscopic structures (Fig.[TjD) 
which are later amplified (Fig. [!}:) by the non-linearities 
in Eq. Furthermore, Fig. [TJi displays a close-up of 
the center of the beam after 7 m of propagation. The pat¬ 
terns are quite similar at both scales. In particular, they 
share the following common features: (i) local maxima 
attracting intensity, depleting the energy around them; 
(ii) strings of intermediate intensity connecting these lo¬ 
cal maxima; (iii) regions of weaker intensity (photon 
bath) around them; (iv) lateral interactions between the 
maxima structures, and (v) the overall shrinking of the 
whole pattern towards a structure with the lower length 
scale. Similar patterns can therefore be observed on two 
spatial scales, two orders of magnitude apart in size. 

Beyond the visual aspect, the self-similarity can be 
quantitatively evidenced by investigating the structure 
factor Sfj) of the laser fluence A = defined by: 


5^(k,r?) = (|i(k,r?)|2), (4) 


z=0m z=7m 



FIG. 1. (Color online) (a), (b), and (c) Evolution of an ini¬ 
tially perturbed fourth order super-gaussian (flat top) laser 
profile of 50 TW at 800 nm, with 5 cm waist. The modu¬ 
lational instability seeds the emergence of a pattern, which 
self-sustains when non-linear effects come into action, (d) 
Magnification of a central zone showing self-similarity. 


increase differs from the monotonic decay of the correla¬ 
tion length that is obtained with thresholded, two-color 
images [l^. At further propagation distances, the flu¬ 
ence clusters either vanish because of dissipation, or get 
squeezed in size because of the energy flux towards their 
center, resulting in a decrease of 


where the hat symbol denotes the Fourier transform, and 
the brackets an ensemble average. Figure [2^ displays 
three spectra corresponding to two stages of the evolu¬ 
tion of the laser beam. Starting from a flat transverse 
spectrum describing the various lengthscales of the initial 
profile modulated with a white noise, the modulational 
instability seeds the emergence of the characteristic pat¬ 
terns at stake here. The peaks on the spectrums after 7 
and 12 m propagation depict the aforementioned multiple 
scales constitutive of the self-similarity. 

Let us define the characteristic length ^ of the trans¬ 
verse patterns using the structure factor by the circular 
average 


^ J S^(k,r))dk 
^ J kS^{k,r])dk' 


(5) 




FIG. 2. (Color online) (a) Structure factor S^{k) calculated 
after 4 m (representative of the initial conditions), 7 m, and 
12 m of propagation. Note that we have suppressed the zero 
peak for clarity reasons, (b) Evolution of the correlation 
length ^ (Eq. ®) for a 800 nm, 50 TW beam of 5 cm waist. 


During the propagation, ^ first increases from the ini¬ 
tial noise correlation length until a maximum length at¬ 
tained at the percolation threshold [l^ (Figurej^b)- This 


The images shown on Figure[T]are reminiscent of many 
models studied by the statistical mechanics community. 
For example, one can note a striking resemblance with 
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coarsening phenomena 4^ ^|. In our case, the trans¬ 
verse low and high intensity regions can be seen as 
two different phases of a generic model evolving under 
Ginzburg-Landau dynamics (48| . 

Such a behavior is generally well reproduced by simple 
spin models with proper time dynamics, and we shall now 
derive lattice spin model (LSM) of filamentation, based 
on the previous observations. The patterns topology we 
aim at reproducing is mainly due to the combined effects 
of modulational instability and Kerr non-linearity. We 
therefore truncate / (Eq. ([2|)) to its cubic contribution 
in Ip. 

The typical patterns shown in Figure [T] strongly sug¬ 
gest to model the electric field by a superposition of nar¬ 
row, Gaussian-like, elementary wavelets. We chose their 
spatial extensions comparable to the lowest-order struc¬ 
ture in the beam, i.e. 10 /rm. Therefore, the field can be 
expanded as 'ip’(r,T]) = A„(r, . 

The universality class unveiled in [1^ suggests that the 
behavior of a lattice model close to criticality should be 
independent from the microscopic detail, and a fortiori 
from the lattice geometry. Hence, we define the set of 
{r„} as a square lattice. Projecting Eq. © on each e''^” 
and identifying real and imaginary parts leads to 


dr,An = -2VAn-V(j)ri-AnA(l)n, ( 6 ) 

AndrjCpn = AAn - AnlVcprif (7) 

“f Aji ^ ^ A^A.fYi cos 0m) • 

£.,m 


By definition, An displays a maximum at r = Tn. Con¬ 
sidering that the phase 0„ is strongly impacted by the B- 
integral 0 , hence by the amplitude An , we assume that 
it also has an extremum at the same location. There¬ 
fore, Eqs. ®-0 can be simplified by cancelling every 
first spatial derivative. As detailed in the Supplementary 
Information [i^ . the spatial self-similarity allows us to 
define the each square lattice site n of area as holding 
two observables, A„ and 0 „, defined as the respective 
averages of the amplitude and phase of the underlying 
small-scale wavelets of the considered cell: 


An — ^ [0n] An^ (8) 

4>n = —2 ^ COS {(pn ~ 4>i) , (9) 

" Wn 

where the dots in (jH])-© refer to a “time” derivative, 
which will reproduce the propagation dynamics of the 
original NLSE, and where n is the discretized Laplacian 
over the four nearest-neighbors. For a site it reads 

Eqs. ©-© are the main result of this Letter. Each 
lattice site can be seen as an individual classical rota¬ 
tor, described by two observables A and cp, which are its 
length and angle, respectively. These rotators evolve un¬ 
der nearest-neighbor interactions, arising from both the 


discretized Laplacians and the last term of Eq. m , ac¬ 
counting for the coarse-grained interference phenomenon. 
For instance, if two lattice neighbors share a common op¬ 
tical phase, their amplitudes will constructively interfere, 
mimicking the situation in which two filaments attract 
each other, and eventually merge. Conversely, if these 
two neighbors feature a relative phase shift of tt, a de¬ 
structive interference will decrease their amplitudes and 
eject their energy towards sites further away, mimicking 
the experimentally observed repulsion [i^ . 

The interaction term Jni = cos(0„ — (pi) is typi¬ 
cal of the spin-glass model, e.g. the soft-spin version 
of the Edwards-Anderson model , character¬ 

ized by the interaction Hamiltonian between spins ai 
H = and evolving under a phenomeno¬ 

logical Langevin equation such as 


d* = -/3 


m 

5a I 


' — P 'y ) Jij aj , 

{j)i 


( 10 ) 


P being the inverse temperature and a Gaussian ran¬ 
dom variable. 

As a test for the relevance of the presented Lattice Spin 
Model (LSM), we will now compare its pattern predic¬ 
tions with the results obtained by integrating the NLSE 
using a standard Split-Step Fourier Method (SSFM). As 
an initial condition, we will use an already slightly prop¬ 
agated beam (by 4 m) with the same properties as in Fig¬ 
ure [T] In general, one could estimate the coarse-graining 
length 5 as being the inverse of the wavelength maxi¬ 
mizing the linear growth rate 7 , since clusters of such 
a size are expected to emerge quicker than others. Do¬ 
ing so with an initial condition as presented here yields 
6 = 924 /rm, in good agreement with the actual physi¬ 
cal range as the size of the photon bath surrounding a 
single filament is typically between 500 and 1000 fim. 
Practically speaking, the coarse-graining leading to the 
definition of the spins (A„, (pn) from the field ip writes as 

|V'(r„-krOpd^r'^ , (11) 

0n = ^ ^ 0(i-n + r')d^r', (12) 

where stands for the lattice cell n of area 5^. Note 
that it is important to average the fluence \ip\'^, and then 
only take the square root instead of directly averaging 
the field ip. This way, Eq. (HU ensures the conservation 
of the photon number Pq = J2n 5^An . 

A direct integration of Eqs. ®-® yields a very good 
qualitative agreement with the reference NLSE patterns 
computed with the SSFM. However, the smallest coarse¬ 
grained length scales of order <5 (typically millimetric) 
behave within the same timescale in the LSM as their 
much smaller counterparts of a hundred micrometers in 
the NLSE hereafter denoted by £c- As a consequence, 
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a pattern arising after a few meters would be predicted 
after only a few centimeters by the LSM. 

However, a linear stability analysis shows that the LSM 
exhibits the same growth rate given by Eq. m as the 
NLSE, which is remarkable. Based on this result, we 
devised a strategy detailed in the supplementary infor¬ 
mation in order to recover the proper dynamics, in¬ 
troducing a rescaled “time” variable r reading 

(13) 

This time renormalization therefore ensures that the 
LSM correctly reproduces the speed of the transverse dif¬ 
fusion of information. 

In our case, we considered a coarse-graining length 
S = 732 /im, i.e. 40 pixels in our reference NLSE nu¬ 
merical resolution. Since we smoothed our initial ran¬ 
dom noise over a length of 4 pixels, we set the small- 
scale cutoff length 4 = d/10 = 73.2 ^m. From Eq. (IT^ . 
we deduce that the time renormalization factor is equal 
to 0.32. Again, this factor lower than unity translates 
the fact that the coarse-graining causes the LSM to act 
on the patterns much quicker than the NLSE does, since 
the former is an upscaled version of the latter. 

We first assess the validity of the LSM by consider¬ 
ing a flat initial phase, namely a real initial ^/>. Figure [3] 
(a-e) compares the fluence pattern obtained from both 
the LSM and the NLSE after approximately 9 m of free 
propagation. Despite the apparent lack of information in 
the initial condition (at 4 m), the LSM remarkably re¬ 
produces the final reference NLSE pattern, showing that 
the interpretation of multiple filamentation in terms of 
interacting spins yields quantitative predictions.This is 
very remarkable as filamentation is generally considered 
as a local phenomenon requiring a high spatial resolution 
in order to capture its salient features. 

To check the fidelity of the phase evolution in a “worst- 
case scenario”, we also considered initial abrupt phase 
jumps from zero to —7r/2 (Figured). This case was cho¬ 
sen because of the difficulty to simulate fields with steep 
gradients, which result in strong diffraction and instabil¬ 
ities that can only be resolved at extremely high reso¬ 
lutions. As a coarse-grained model intrinsically cannot 
capture such features, such a situation should check the 
robustness of our lattice spin description. 

We modulated the aforementioned amplitude pattern 
by a phase mask displaying the word “unige”. Figure [3] 
(g), (h), and (j) shows a remarkable agreement after 3 m 
of free propagation, despite a glitch on the reproduction 
of the letter “g”. 

Moreover, these results are not tributary to a fine- 
tuned choice of the coarse-graining length 6. One can 
freely choose it in the aforementioned physically accept¬ 
able range and still obtain reasonable results, whereas a 
decrease of resolution in the SSFM method rapidly leads 
to erroneous simulations. 
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FIG. 3. (Color online) Comparison between the lattice spin 
model and the coarse-grained result of the NLSE simulation 
using the SSFM. (a-e), fiat initial phase; (f-j), worst-case sce¬ 
nario with n/2 phase jumps, (a), (f) Initial condition for the 
LSM, originating from a full resolution SSFM integration; (b), 
(g) LSM output; (c), (h) NLSE output using SSFM; (e), (j) 
Horizontal cuts across the model ouptuts. 


These two test cases highlight the relevance and even 
the quantitative accuracy of the LSM. For smooth ini¬ 
tial phases, the relative error on intensity stays below 
10%. The importance of nearest-neighbor interaction was 
further demonstrated by switching off the corresponding 
term in Equation (|9]). The beam then keeps a smooth 
shape very different from the self-structuring of the beam 
observed in both experiments and NLSE simulations. 

It is quite straightforward to derive richer lattice mod¬ 
els encompassing more phenomena, such as e.g. other 
non-linearities, saturation mechanisms, plasma genera¬ 
tion (see Eq. @) or even air turbulence. This would 
simply require to expand the additional physical model 
on the wavelet basis, and then simplify all the remaining 
terms by keeping in mind the nearest-neighbor picture. 

As a conclusion, we took advantage of the self¬ 
similarity of multiple filamentation patterns to introduce 
the description of laser multiple filamentation as a Lat¬ 
tice Spin Model with glassy-like dynamics. The numeri¬ 
cal benchmarks showed an excellent agreement with the 
full calculations, demonstrating the robustness of such 
a novel interpretation, that can also be related to the 
recent observation of a percolation-like phase transition 
in such a system [l^ . Furthermore, as a consequence 
of the coarse-grained description, the small lattice sizes 
at play allow computing times faster by two orders of 
magnitude as compared to standard SSFM calculations. 
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Such a speed-up opens the way to statistical studies of, 
e.g., beam propagation through turbulence, or explicit 
inversion of non-linear Lidar measurements 53“52| of at¬ 
mospheric trace constituents. 

In a wider scope, our approach only relies on the struc¬ 
ture of the NLSE, not on a particular nonlinearity (i.e. 
a particular function /), nor its application to a specific 
physical system. Therefore, it can be generalized to other 
fields of physics described by the NLSE, where such self¬ 
similarity could also be observed and exploited ^55|. 

We wish to warmly thank an anonymous referee for 
very valuable comments and suggestions. We grate¬ 
fully acknowledge fruitful discussions with T. Giamarchi, 
M. Brunetti and S. Hermelin. We acknowledge financial 
support from the European Research Council Advanced 
Grant “Eilatmo” and the Swiss National Science Eounda- 
tion (Grant 200021-155970). 


DIMENSIONLESS UNITS 

Considering the non-linearity function/(I'0p) = IV'P — 
where K is the number of photons 
required for medium ionization, the dimensionless units 
read 

= (14) 

-tp = (15) 

X = {'y/{akon2)^) (16) 

where x (transverse coordinate) and z (propagation dis¬ 
tance) are in meters, A (the electric field envelope) in 
V.m“^. fco = 27r/Ao is the central wavenumber in m“^, 
for which we use Aq = 800 nm. The non-linear refractive 
index n 2 is set to 1.2 x 10“^^ m^.W“^. The parameters 
a and 7 describe the medium delayed response, and can 
be found in reference (s^ . 


DERIVATION OF THE LATTICE MODEL 
EQUATIONS 

In this section, we shall explicit the derivation of the 
lattice model starting from the observation that the elec¬ 
tric held Ip can be developed on a basis of Gaussian 
wavelets, as emphasized in Figure 0] Practically speak¬ 
ing, we expand the held as 

p^ir,rj) = (17) 

n 

Injecting this decomposition into the NLSE, one obtains 
the following system: 

dr,An = - 2 VA„ • V(pn - A„A^„, (18) 

Andr,(pn = l^An - An\V (pn\^ (19) 

T Arf^ ^ ( A^Ayji cos (^(pi (pm) ■ 

.^,771 



ri 12 13 14 Is re 17 rg rg r^o 


FIG. 4. Principle of the wavelet decomposition and lattice 
discretization procedure. Each wavelet k of the decomposition 
is supposed to be centered on rk. 


The extrema of A„ and (pn at each beam center r = rn 
allow to cancel out every hrst spatial derivative at these 
locations. The system (IT^ - (l2nil then rewrites as: 

5^A„ = -(A<^„)A„, (20) 

drjcpn = ^ A^(r„)Am(r„) (21) 

" e,m 

X cos [(pe{rn) - (/'m(rn)] ■ 

As the width of the individual wavelets has been chosen 
comparable with that of the lattice cells, interactions be¬ 
tween elementary beams can be neglected except for the 
four nearest neighbors, greatly simplifying the last term 
in Eq. (1^ : 

A A 

dr,(Pn = ^+Al ( 22 ) 

+An ^ A^(r„) cos [(pn - (piirn)], 


where the notation {i)n denotes a summation over the 
nearest-neighbors on the square lattice of the elemen¬ 
tary beam. Note that the second term in the right-hand 
side of Eq. (|22|l corresponds to the standard 712 / cumula¬ 
tive phase shift term, also known as the R-integral. This 
phase evolution equation is similar to a short-range dis¬ 
ordered Kuramoto model a paradigmatic model 

for synchronization. 

We shall now take benefit of the self-similarity of our 
problem to upscale the model from the 100 /im micro¬ 
scopic filamentary structures, to that of the millimeter 
sized aggregates. In that purpose, we define a new lat¬ 
tice, with coarser cells of width hereafter denoted by S. 
To get rid of the continuous spatial representation r in 
Eqs. (|^ and (l22)l . we replace the Laplacian operator 
A by its square-lattice discretized counterpart, named n, 
which action on a site n = {i,j) reads 


^ + (l^i-i,j + 4(1,j+i + (pi, 3 - 1 ) ■ 

(23) 

We also replace the wavelets defined by ^^(rn) (resp. 
</>f(rn)) by A^(rf) (resp. (pi{vi). Note that this is a strong 
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assumption, but essential for the final nearest-neighbor 
interacting model, which finally reads 

A.n = K [(/>n] 

(j)n = —2^ ^ ^ COS {(j)n ~ 4>i) ■ (25) 

-^n ,,, 


TIME RENORMALIZATION 


In this section, we shall demonstrate the correct time 
rescaling so as to obtain the correct propagation dynam¬ 
ics with the lattice model. Neglecting the last nearest- 
neighbour coupling in equation (jH]), we find that like for 
the NLSE, the plane-wave defined by the homogeneous 
amplitude A = A* and phase (j)* = (j) + A*^t is a steady- 
state solution. Furthermore, linearizing the system (IMll - 
(051) and considering the first-order corrections of 6A and 
54> yields, in the Fourier space. 


5A 

dt \ 5^ 


0 ^ 

2A* - ^ Q )\5(t>) 


(26) 


The eigenvalues of the Jacobian matrix are given by 

7 = ±kV2A*^ - P, (27) 


In Fourier space, the spatial averaging acts as a low- 
pass filter, and brutally cuts all spatial frequencies above 
kc = 21115. As a result, the growth rates calculated for 
both the original pattern and its coarse-grained counter¬ 
part might differ because of the change in the spectrum 
distribution, leading to an erroneous description of the 
dynamics by the lattice model. In order to tackle this 
issue, we rescaled the time variable in the lattice model 
so as to recover the same apparent growth as for the full 
equation. 

More specifically, we rescaled time in the coarse¬ 
grained case so that the Fourier-averaged initial time 
derivatives of the growth rate coincide for both cases. 
Namely, we look for a time variable r satisfying 



ip{k, r])dk 


rj—O 


A 

dr 


j A{k,T)dk 


(28) 


This approach is motivated by the fact that we want the 
lattice model to reproduce the dynamics of the NLSF al¬ 
ready in the beginning of the propagation. For instance, 
we chose to rely more on the time derivative at the origin 
rather than the maximum growth rate. 

Plugging the fields’ expressions in the linear regime in 
Fq. (l28l) . we obtain 


^{k,r] = 0)j{k)dk = / A{k,T = 0)^{k)dk, (29) 


In order to express r with respect to 77 , let us consider a 
full resolution initial condition ipQ exhibiting a flat spec¬ 
trum, with a high-frequency cutoff 27r/£(,. In this case, 
the structure factor reads = <S'oX[o, 27 r/ 4 ]) where xe is 
the indicatrix of the set E. This corresponds to a plane 
wave perturbed by a white noise with a correlation length 
of Ic- In virtue of Parseval’s theorem, the spectral power 
amplitude relates to the initial power Pq by 

^0 = — ^ 0 , (30) 

After the coarse-graining procedure, still in virtue of Par¬ 
seval’s theorem, we have 


Y,5^\Ak\^ = Po, (31) 


where the summation is performed over the reciprocal 
lattice, A being the Fourier-transformed coarse-grained 
field. Since the frequencies higher than 2it/5 have been 
wiped out from the new spectrum by the coarse-graining, 
Fq. pil) implies that the new spectral power ampli¬ 
tude must read 

S^o^ = Soj-. (32) 

Since 6 > £c, the new spectrum displays a higher ampli¬ 
tude in order to ensure the photon number conservation, 
despite the disappearance of high-frequency modes, that 
corresponded to the microscopic detail that we got rid of 
by spatial averaging. 

It is now our aim to calculate the average growth 
rate ( 7 ) for both cases, given the two different spectral 
probability distributions 

fik) = y^^X[o, 2,7/4], (33) 

2 , 7 / 5 ] ■ (34) 

Let us calculate the probability distribution of the re¬ 

striction of 7 to its real positives values. The probability 
of observing a growth rate lower than 7 for a random k 
chosen with probability / (which can be given either by 
Fqs. (|33l) or (|34|) 1 then reads 


pAV2 

P{y < 7) = / /(fc)dfc+ / ^ f{k)dk. 

Jo J 

(35) 

After differentiation with respect to 7 , one obtains the 
probability distribution function for the positive growth 
rate. It is then straightforward to understand that the 
coarse-grained average growth rate relates to its original 
counterpart through 


where the growth rate 7 is given in each model by 

Fq. (EH). 


(^cg) 



(36) 
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which gives the final relationship between the time vari¬ 
able r in the lattice model and the original distance vari¬ 
able rj in the NLSE, 



COMPUTATIONAL CONSIDERATIONS 

Given a N x N square lattice, the model integration 
is of complexity 0{N^), while the SSFM, helped by the 
Fast-Fourier-Transform algorithms available at hand, has 
a complexity scaling as 0{N^\ogN), which is slightly 
greater. 

The result for 9 meters of propagation for a realistic 
field as presented in the main article is obtained after 
only 35 s of computational time on a desktop computer, 
while the NLSE took 5 h on a dedicated workstation, 
without GPU acceleration. 

Even if the recent GPU accelerators can greatly im¬ 
prove the computational performance of the SSFM, the 
resolution and complexity of the lattice model are far 
too low to compare. The lattice model provides a new 
paradigm for laser multiple filamentation, the computa¬ 
tional speed-up is a consequence of such a simplification, 
though at the expense of the microscopic detail. 
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